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Abstract 

This paper describes the numerical simulation of the NASA Langley Research Center supersonic H2-M1: 
combustion chamber performed using two approaches to model the presumed probability density function 
(PDF) in the flamelet progress variable (FPV) framework. The first one is a standard FPV model, built 
presuming the functional shape of the PDFs of the mixture fraction, Z, and of the progress parameter, A. In 
order to enhance the prediction capabilities of such a model in high-speed reacting flows, a second approach 
is proposed employing the statistically most likely distribution (SMLD) techcnique to presume the joint 
PDF of Z and A, without any assumption about their behaviour. The standard and FPV-SMLD models 
have been developed using the low Mach number assumption. In both cases, the temperature is evaluated 
by solving the total-energy conservation equation, providing a more suitable approach for the simulation 
of supersonic combustion. By comparison with experimental data, the proposed SMLD model is shown to 
provide a clear improvement with respect to the standard FPV model, especially in the auto-ignition and 
stabilization regions of the flame. 

Keywords: Joint Presumed PDF modelling, Hydrogen-Air combustion, Reynolds-Aver aged Navier-Stokes 
equations 


N omenclat ure 

C Progress variable 
D(f) Diffusivity of 0 
H Total enthalpy 
k Turbulent kinetic energy 
Ma Mach number of the flow 
p Pressure 

P{x) Density-weighted probability density func¬ 
tion 


PsML,2{x) Density-weighted probability density 
function evaluated as the statistical most likely 
distribution with the second order moments 

q Heat flux 

Qreact Heat release rate 
R Gas constant 

Re Reynolds number of the flow 
T Temperature 
Tu Turbulence intensity 
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u Velocity 

Mass fraction of (p 
Z Mixture fraction 
/3{x) /^-distribution 
r Euler function 
6{x) Dirac distribution 
A Progress parameter 
/i Dynamic viscosity 
jXx Lagrangian multiplier 
u Kinematic viscosity 
p Density 

(j) Generic thermo-chemical quantity 


(j) Favre-averaged value of (j) 

(f)" Fluctuation of 0 in the Favre-averaging process 

0"^ Variance of (j) 

(j) Reynolds-averaged value of (j) 

(p' Fluctuation of (j) in the Reynolds-averaging pro¬ 
cess 

^ Error function 

X Scalar dissipation rate 

Xst Value of y at the stoichiometric point 

uj Turbulent kinetic energy specific dissipation rate 

uj(j) Production therm of (j) 


1. Introduction 

In the last years, the development of propulsion systems based on air-breathing engines for super-/hyper- 
sonic flying vehicles has fostered the study of supersonic combustion, see, e.g., the review of Cecere et ah [T]. 
In these systems, hydrogen is one of the preferred fuel because of its properties in terms of very short ignition 
delay time and high energy per unit weight. The investigation of hydrogen supersonic combustion presents 
significant difficulties and high costs either following the experimental approach or the numerical one. In 
fact, in supersonic combustion, the mixing time scales are comparable to i^ 2 -Air reaction time scales [2]. 
Moreover, high-Reynolds-number turbulent combustion is a formidable multi-scale problem, where the in¬ 
teraction between chemical kinetics, molecular, and turbulent transport occurs over a wide range of length 
and time scales. These features pose severe difficulties in the analysis and comprehension of the basic phe¬ 
nomena involved in supersonic combustion. 

Concerning the numerical approach, in recent years, the need for efficient tools having affordable compu¬ 
tational costs has driven the research towards: i) studying turbulent combustion in order to understand 
the interaction between turbulence and chemistry mUEli]; ii) improving kinetic schemes to describe the 
combustion process 0019 ]. Higher accuracy can be achieved by employing models based on detailed kinetic 
mechanisms, but this usually leads to prohibitively expensive calculations. Therefore, reduced models are 
often employed to condensate the reaction mechanisms and cut down the computational costs [7] . Simplified 
approaches to combustion modelling have been proposed to further reduce the number of equations to be 
solved; for instance, the reduction of the chemical scheme in intrinsic low dimensional manifolds (ILDM) [10]; 
the flamelet-based approaches such as the fiamelet-progress variable (FPV) [H] or flame prolongation of 
ILDM (FPI) [12]; and fiamelet generated manifolds approach (FGM) [T3] . 

The present work is based on the FPV model for non-premixed flames. Standard steady FPV models 
are built under the low Mach number hypothesis [Hiigiiaiiii, computing the combustion process at a 
uniform given pressure. Obviously, in supersonic combustion, density variations due to the dynamics of 
the flow cannot be neglected. For this reason, here we employ the modified FPV approach proposed by 
Oevermann [18], and recently employed by other authors [19], where the fluid temperature is not obtained 
from the fiamelet libraries as in the standard model, but is calculated by solving the full set of conservation 
equations for compressible flows, including the energy equation. More recent works have proposed an 
extension of FPV combustion models to compressible flows by using a source term linearly rescaled with the 
pressure value and a perturbation of the low-Mach number fiamelet accounting for compressibility [iniini. 
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For the case of non-premixed combustion of interest here, mixing must bring reactants into the reaction 
zone so as to activate and maintain the combustion process. Such flames are characterized by a local balance 
between diffusion and reaction [22]. The basic assumption of the flamelet model is that each element of 
the flame front can be described as a small laminar flame, also called flamelet. Therefore, for a steady 
flow, the flame structure can be described as a function of the mixture fraction, Z, and of the progress 
parameter, A m- For turbulent combustion, a probability density function (PDF) is needed to compute 
the mean value and the variance of the thermo-chemical variables. The definition of such a PDF is critical 
since it has a strong impact on the solution. The aim of this work is to study the applicability of the 
statistically most likely distribution (SMLD) [24] approach to model joint-PDF of Z and A in the case of 
supersonic combustion. The proposed joint-SMLD approach is very interesting since it represents a good 
compromise between computational costs and accuracy level. The results obtained using the proposed model 
are validated versus experimental data and compared with numerical results available in the literature as 
well as with those obtained using the standard FPV model. 

This work is organised as follow: sections and provide the theoretical description and some details 
of the numerical discretization for the two models. The comparison between their numerical results for the 
simulation of the NASA Langley Research Center supersonic hydrogen flame are presented in section]^ along 
with reference numerical and experimental data. Finally, some conclusions are provided. 


2. The flamelet progress variable models 


For the case of non-premixed combustion of interest here, the basic assumptions of the flamelet model are 
fulfilled for sufficiently large Damkohler number. Da. In fact, when the reaction zone thickness is very thin 
with respect to the Kolmogorov length scale, turbulent structures are unable to penetrate into the reaction 
zone and cannot destroy the laminar flame structure. Effects of turbulence only result in a deformation 
and straining of the flame sheet and locally the flame structure can be described as function of the mixture 
fraction, Z, the scalar dissipation rate, y, and the time. The scalar dissipation rate, y = 2Dz(VZfl, is a 
measure of the gradient of the mixture fraction representing the molecular diffusion of the species in the 
flame, Dz being the molecular diffusion coefficient of the chemical species. Therefore, the entire flame 
behaviour can be obtained as a combination of solutions of the laminar flamelet equation. In the present 
work we consider a further simplification assuming a steady flamelet behaviour, so that chemical effects are 
entirely determined by the value of Z, whereas y describes the effects of the flow on the flame structure 
according to the following steady laminar flamelet equation (SLFE) for the generic variable 0: 


X d'^4> _ . 
^2dZ‘^ 


( 1 ) 


In equation Q, p is the density and is the source term related to (j) [53], different from zero in the case 
of finite rate chemistry. In particular, in this work, the EPV model proposed by Pierce and Moin mills] is 
employed to evaluate all of the thermo-chemical quantities involved in the combustion process. This approach 
is based on the parametrization of the generic thermo-chemical quantity, 0, in terms of the mixture fraction, 
Z, and of the progress parameter. A, instead of y: 


cP = F^{Z,A). (2) 

Using such a parameter, independent of the mixture fraction, one can uniquely identify each flame state 
along the stable and unstable branches of the S-shaped curve. A suitable definition of A leads to a dramatic 
simplification of the presumed PDE closure model. On the other hand, the solution of the transport equation 
for A is quite complex since it requires non-trivial modelling of several unclosed terms [25]. In order to 
overcome such a difficulty, the progress parameter is derived from a reaction progress variable, (7, such as 
the temperature or a linear combination of the main reaction products, whose behaviour is governed by 
a simpler transport equation. Therefore, a transport equation for C is solved and the flamelet library is 
parametrized in terms of Z and C. Requiring that the transformation between A and C be bijective, from 
equation 0 one has 

A = Fc\Z,C), 
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(3) 



and any thermo-chemical variable can be expressed as: 

^ = F4,{Z,Fc\Z,C)). (4) 

The choice of the progress variable is not unique and some recent works discuss in details this issue proposing 
a procedure for its optimal selection [26l|271|28]. A suitable definition for the progress variable is the sum 
of the mass fraction of the main products [26]; for hydrogen combustion m- 

c = Yh,o- (5) 


A stretching with respect to the minimum and maximum conditioned value of C over the mixture fraction 
is made, leading to the following form of the progress parameter: 

C — CMin\Z 


A = 


( 6 ) 


CMax\^ CMin\^ 

Equation ([^ is taken as the solution of the SLFE Q . Even if cases in which there is not a unique mapping 
of this solution in function of Z and A are reported in the literature mills], such events are excluded from 
the solution family composing the flamelet library, since they are very close to the equilibrium limit. Since 
flamelet libraries are computed in advance and are assumed to be independent of the flow field, one has to 
model the dependence of y on Z [ISl[29l|30]. In this work, the functional form of x(^) has been taken from 
an idealized flow configuration, as proposed by Peters ED; the distribution of the scalar dissipation rate in 
a counterflow diffusion flame, ^(^), is employed, scaled in the following way: 


x{Z) = Xst- 


(7) 


where Xst and Zgt are evaluated at the stoichiometric point m- 

Erom equation Q one can derive the generic thermo-chemical quantity (j). When a turbulence model is 
used, the Eavre-averaged value of (j) and its variance are given as: 


(t) = J J F^Z, A)P{Z, A)dZdA, 

^ = J j iP^Z, A) - ^YP{Z, A)dZdA. 


In the above equations, P{Z^ A) is the density-weighted PDE, 

P 

P{Z, A) is the joint PDF and p is the Reynolds-averaged density. As usual, 0 can be decomposed as: 


0 0 + < 


P 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 


and 


P — pPp'i (12) 

where 0" and p' are the fluctuations. The joint-PDF, P{Z,A), plays a crucial role in the definition of the 
model, affecting both its accuracy and computational costs. Moreover, the choice of such a function is 
not straightforward because of the unknown statistical behaviour of the two variables Z and A [25] and its 
definition is still an open problem whose solution is being pursued by several researches [32] |33j |34l ES] . 
The aim of this work is to validate a more general model based on the statistically most likely distribution 
(SMLD) [24] for the joint PDF of Z and A [36]. The performance of such a combustion model are assessed by 
computing a hydrogen-air supersonic flame and comparing the results with those obtained using the standard 
FPV model. In the present work, the Reynolds-Averaged Navier-Stokes equations with k-uj turbulence 
model are solved and both the standard-FPV and the FPV-SMLD models employ the total energy 
conservation equation to evaluate the temperature field in order to improve the simulation of compressible 
reacting flows. 
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2.1. Presumed probability density funetion model 

In this section, the standard FPV model [23] (called here model A) and the FPV-SMLD model (called 
here model B) [36| are briefly described. 

The evaluation of the average quantities in equations and © requires the PDF to be known or 
somehow presumed. Such a PDF establishes the statistical correlation between Z and A. Employing the 
Bayes’ theorem, 

P{Z,K)=P{Z)P{k\Z), (13) 

one usually presumes the functional shape of the marginal PDF of Z and of the conditional PDF of A\Z. 
In model A, the basic assumption is the statistical independence between Z and A, so that P(Z, A) = 
P{Z)P{A). Furthermore, the statistical behaviour of the mixture fraction is described by a ^^-distribution. 
In fact, even though the definition of P{Z) is still an open question [24|, it has been shown by several authors 
that the mixture fraction behaves like a passive scalar whose statistical distribution can be approximated 
by a /d function [38l |39l [40] . The two-parameter family of the /^-distribution in the interval x G [0,1] is given 
by: 

/3(a;;x,a;"2) =a;“-i(l(14) 

where T{x) is the Euler function and a and b are two parameters related to x and x"‘^ 


c{x — x^ — x"‘^) 

rf.n2 


b = 


(1 — x)(x — x‘^ — x"‘^) 


p//2 


(15) 


Moreover, P(A) is chosen as a Dirac distribution, implying a great simplification in the theoretical framework. 
With these assumptions, the Eavre-average of a generic thermo-chemical quantity is given by: 


(t) = J j F4,{Z, K)j3{Z)5{k-K)dZdC = J F^Z, K)j3{Z)dZ. 


(16) 


Therefore, in addition to the conservation and turbulence model equations, one has to solve only three 
transport equations (for Z, Z"‘^ and C) to evaluate all of the thermo-chemical quantities, thus avoiding the 
expensive solution of one transport equation for each chemical species. 

Model Debased on the SMLD approach to model the joint PDE, does not need any assumption about 
the form of P{Z^ A). Eollowing such an approach, the probability distribution can be evaluated as a function 
of an arbitrary number of moments of Z and A. It is noteworthy that, even though equation 0 is based on 
the assumption that Z and A are independent, one can properly take into account the statistical correlation 
between Z and A employing the SMLD joint-PDE in the evaluation of the effects of turbulence [32] . 

In this work, the first two moments of the joint probability density function P(x), where x = {Z^K^^ 
are assumed to be known; therefore, the joint-PDE reads [36] : 

PsmlA^^^) = ^exp|-Li,i(Z-Z)+/ii,2(A-A) 

k2M{Z - Zf + M2,12(^ - Z){k - A) + /X2.2i(A - A)(Z -Z)P M2.22(A - A)^] }. (17) 

In the equation above, /io is a scalar, /ii is a two - component vector, and is a square matrix of rank two: 


1 r 

2 


/io 

^kl ~ 1^2,kn 


= J dxPsML,2{x), 

= J dxd,XsML,2{x) = ^(1; I, ) - ^(0; I, ef), 

= J dxdxk{{xi - ii)PsML,2{x)) = ld{l\ik,i'kei) - 
5 


(18) 

(19) 

( 20 ) 








where and I indicate the vector components; and are the mean {xi)^ the fluctuation {xi—Xi) 

and the variance {x'p) of the i-th component of x, respectively; finally, P indicates the beta distribution 
function. In model B, one has to solve four additional transport equations (for Z, C, and C"‘^) to 

evaluate all of the thermo-chemical quantities. 


3. Governing equations 

3.1. Flow equations and numerical solution procedure 

The numerical method developed in [37] has been employed to solve the steady-state RANS equations 
with k-uj turbulence closure. For an axisymmetric multi-component reacting compressible flow the system 
of governing equations can be written as: 

dtQ + d,{E - E,) + dy{F - F,) = S, (21) 

where t is the time variable; x and y are the axial and the radial coordinate, respectively; Q={p, pu^^pUy^pH— 
Pt^~pk^ puj^pRn) is the vector of the conserved variables; E^ F, and Ey^ Fy are the inviscid and viscous flux 
vectors [41], respectively; S is the vector of the source terms; p, {ux^Uy), H indicate the Reynolds-averaged 
yalue^of density, the Favre-ayeraged values of velocity components and specific total enthalpy given 
H = h-\- ^{u‘lEUy) E^k with h accounting for the species enthalpy per unit mass, respectively; pt = p + |/c, 

p being the Favre-averaged value of pressure; k and uj are the Favre-averaged values of the turbulence kinetic 
energy and of its specific dissipation rate; Rn is a generic set of conserved variables related to the combustion 
model. In this framework, R^ is the set of independent variables of the flamelet model, namely, Z, Z"^, C, 
(see the following subsection). 

The heat flux in the total energy equation, namely. 


Ns 

q = pCpDj^S/T T ^ ^ pV^iYnh-n E qreacti 
n=l 

is composed of three terms since the Dufour effect is neglected; Dt is the thermal diffusivity, Cp the specific 
heat at constant pressure; the mass diffusion term is modelled by the Fick law considering Vn = — 
assigning mixture diffusivity, F, to each species and so assuming unitary Lewis number [22]. The third term 
at the right hand side of the above equation represents the heat release rate: 

N 

Qreact — ^ ^ 
k=l 

where k is the species index, Ahj ^ is the mass formation enthalpy and Ek is the production rate of species 

k. 

3.2. Turbulent FPV transport equations 

For the case of turbulent flames, the solution of the SLFE, namely equation is expressed in terms of 
the Favre averages of Z and C and of their variance. Using model A, one can tabulate all chemical quantities 
in terms of Z, Z"‘^ and (7, since, due to the properties of the (^-distribution, the model is independent of 
C"‘^. On the other hand, model B expresses (j) also in terms of C"‘^ and therefore an additional transport 
equation needs to be solved. In this case, the transport equations for the combustion model (included in 
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equation (HI)) are written as: 


dt{pZ) + V • (puZ) 
dtipC) + V • (p^c) 


V- 

{D + D*~)pVZ 

5 


(24) 

V- 

\d + d^~)pVZ"^' 

-pX + 2pD|(VZ)^ 

(25) 

V- 

{D + Dy)pVC 

+ pwc, 

(26) 

V- 

'{D + D^~)pVC''^' 

- pxc + 2pDy{VCf + 2pOtPl;, 

(27) 


where xc is modelled in terms of Z"‘^ and C"‘^ 


namely xc 


D 


vjPr is the diffusion 


coefficient for all of the species; v and Pr are the kinematic viscosity and the Prandtl number, respectively; 
D^~ = = J^tlSct are the turbulent mass diffusion coefficients, Sct being the turbulent 

Schmidt number equal to 0.8; finally, ujc is the source term for the progress variable precomputed and 
tabulated in the fiamelet library. At every iteration, the values of the fiamelet variables are updated using 
equations (24)-(27) and the Favre-averaged thermo-chemical quantities are computed, using equation 
Such solutions provide the mean mass fractions which are used to evaluate all of the transport properties of 
the fluid, namely the molecular viscosity, the thermal conductivity and the species diffusion coefficients. 
The evaluation of the fiamelet library has been performed using the detailed kinetic scheme proposed by 
Saxena and Williams m- 244 sub-reactions upon 50 species. The fiamelet library has been computed over 
a grid with 250 uniformly distributed points in the Z and C directions, and 50 uniformly distributed points 
in the Z"‘^ and directions and a quadri-linear interpolation scheme is used. The fiamelet library has 
been evaluated considering a constant background pressure equal to 100 kPa and boundary conditions for 
Z a^d C corre^onding to the conditions of the air and H 2 strums. Indeed, the air stream is represented 
as Z = 0 and (7 = 0, while the fuel jet is given by Z = 1 and (7 = 0. The solution of the SLFE has been 
obtained by using the FlameMaster code [43] . 


3.3. Numerical scheme and boundary conditions 

A cell-centred finite volume space discretization is used on a multi-block structured mesh. The convective 
and viscous terms are discretized by the third-order-accurate Steger and Warming [44] flux-vector-splitting 
scheme and by second-order-accurate central differences, respectively. An implicit time marching procedure 
is used with a factorization based on the diagonalization procedure of Pulliamm and Chaussee [45], employing 
a scalar alternating direction implicit (ADI) solution procedure [46]. Steady flows are considered and the 
ADI scheme is iterated in the pseudo-time until a residual drop of at least five orders of magnitude for all of 
the conservation-law equations (21) is achieved. Characteristic boundary conditions for the flow variables 
are imposed at inflow and outflow points. In particular, a plug flow is imposed at the inlet points of the 
computational domain so as to match the experimental conditions at the inlet section of the chamber, see 
table Moreover, /c, cj, and Rn are assigned at inflow points, whereas they are linearly extrapolated at 
outflow points. Finally, no slip and adiabatic conditions are imposed at walls, where k is set to zero and uj 
is evaluated as proposed in [47] : 


cj = 60 


0 . 09^2 


(28) 


n,l 


where is the distance of the first cell center from the wall; the homogeneous Neumann boundary condition 
is used for Rn (non-catalytic wall). Finally, symmetry conditions are imposed at the axis. 


4. Results and discussion 
4 . 1 . Description of the test case 

The hydrogen-air supersonic combustion burner studied by Cheng et al. [48| at the NASA Langley 
Research Center has been considered as a suitable test case for the proposed method. At the inlet section 
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Figure 1: Schematic of the Cheng’s Flame burner. The red-dashed line represents the exit section of the burner and so the 
inlet section of the chamber. 


of the chamber, the supersonic flow (see figure is characterized by an annular axisymmetric jet of hot 
vitiated wet air at Mach number equal to 2, average axial velocity of 1420 m/s, temperature of 1250 K 
and pressure of 107 kPa. The values of the inner and outer diameters of the vitiated air annular duct 
are equal to 3.812 mm and 17.78 mm, respectively. The composition of the vitiated air is: Yq^ = 0.201, 
YH 2 O = 0.255 and = 0.544, produced by a pre-combustion at low temperature [48]. The diameter of the 
fuel stream is d^ef = 2.362 mm, taken as the reference length. The hydrogen flow is estimated to be chocked 
with average axial velocity of 1780 m/s, temperature of 545 and pressure of 112 kPa. The operating 
conditions reported by Cheng et al. [48] are summarized in table together with the turbulent intensity 
levels, Tuin^ and the turbulent viscosity, The computational domain, shown in figurej^ is axisymmetric 
and includes the divergent part of the air nozzle, necessary to recover the correct flow quantities at the inlet 
section of the chamber; it extends 150 dref and 50 dref along the axial and radial directions, respectively, 
and, after a grid refinement study, has been discretized using about 100000 cells. The characteristic cell 
lengths at the exit of the divergent part of the nozzle are about 0.13 mm and 0.05 mm in the axial and 
radial directions, respectively. Numerical results obtained using the two models described in the previous 
sections are discussed and validated versus the experimental data of Cheng et al. [48] and the numerical 
results of [11149]. In particular, a detailed analysis of the flame structure is provided. 

4 . 2 . Comparison between numerieal and experimental data 


Table 1: Parameters for the simulation of the Cheng’s combustion chamber m 


H 2 jet 

Ma 

1 

u{m/s) 

1780 

T(K) 

545 

p (kPa) 
112 

Yh 2 

1 

To. 

0 

YH 2 O 

0 

Yn 2 

0 

5.78% 

0.00023 

Air stream 

2 

1420 

1250 

107 

0 

0.201 

0.255 

0.544 

10.24% 

0.00034 


Figure [^provides the Mach number contours obtained using model A and model B. The close-up view 
of the near-burner region (bottom-left and -right panels) indicates that the inlet Mach number value is very 
close to the experimental data for both computations. The distributions of the streamwise component of 
the velocity, iz, at several abscissae along the chamber are shown in figure [^ A good agreement with the 
experimental data [48] is obtained for both models. Figure (top frame) shows a qualitative comparison 
between the temperature contours evaluated with the two models. It appears that the computed flame 
shapes are quite different. Model A predicts a reaction zone attached to the burner, whereas model B 
correctly predicts the flame detachment. In the bottom frame of figure [^ one can find the two corresponding 
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Figure 2: Computational domain discretization (lower frame) and a detail of the injectors (upper frame). 



(a) Model A (b) Model B 


Figure 3: Mach number contours obtained with model A (left) and model B (right), with a close up of the near-burner region 
(bottom). 
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Figure 4: Stream-wise velocity distributions at several radial sections: model A, red dashed line; model B, black solid line; 
symbols, experimental data. 
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progress-variable contours. The solution obtained using model A provides non-zero values for the progress 
variable in the region close to the burner, indicating that the reaction is active. Figure presents the 
temperature field obtained using model A (left panel) and model B (right panel) in comparison with the 
results of an accurate LES by Boivin et al. [7]. One can observe that model B predicts the lift-off height, 
evaluated as the position of maximum temperature gradient, at about x = 26 dref^ whereas model A predicts 
an height of about 9 d^ef- The results by Boivin et al. [7], used as reference, predicts the stabilization at 
25 dref- For both models, the jet shape and the aperture of the flame at x/d^ef — 50, equal about to 10 
are in good agreement with the experimental data (not shown). For a quantitative analysis of the results, 
figures [7}|^ provide the comparison among the results obtained using the two FPV models, the numerical 
results provided in m sg, and the experimental data |48]. Figure shows the distribution of the main 
thermo-chemical quantities along the axis of the burner. The evaluation of the Z and are very slightly 
improved by model B; a reasonable improvement of the prediction of the distribution of OH mass fraction is 
also obtained; whereas, H 2 and H 2 O mass fractions are quite well predicted in comparison with model A. It 
appears that the flame core is not well reproduced using model A, which indeed provides a too high reaction 
rate, so that the reactions may occur wherever the two flows (i72“Wet-air) are mixing. As a consequence, 
model A predicts Yh^o maximum at about 10 x/dref, as shown in the middle right frame of figure]^ On 
the other hand, model B can reproduce the ignition spatial delay with a strong increase of the temperature 
at about xjdref = 30 (figuremiddle-left frame). 

Figure also provides the numerical results obtained with an LES by Boivin et al. [7], for comparison; 
they appear to be in good agreement with the results obtained by model B, confirming the good level 
of accuracy achieved by such a model. Analysing the temperature and water mass fraction experimental 
distributions shown in figure]^ one can notice that the temperature begins to decrease at x/dref = 50 
whereas the water mass fraction increases up to about xjdref = 65. This is due to the entrainment of 
the surrounding unburnt-air into the flame, which induces an increase of the mass flow rate and allows the 
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Figure 6: Comparison between the temperature contours computed by model A (left), Boivin et al. [?] (center), model B (right). 
The contour lines are plotted from lOOOiF to 2500iF with step equal to 250iF. 


reactions to occur even if the flame cannot heat up the surrounding flow. In particular, the early increase 
of the temperature and of the water mass fraction, obtained by model A, indicate that the heat is released 
very close to the burner. As a consequence, the entrainment of the surrounding unburnt-air into the flame 
is strongly overestimated using such a model with respect to model B. The entrainment can be evaluated 
by computing the total mass flow rate through the inlet of the chamber, at xjdref = 0, and at the sections 
xjdref = 25 and xjdref = 50. Due to the axial-symmetry, the effective surfaces intersected by the ffow 
can be estimated as the circles with radius equal to the maximum extension of the temperature profiles, 
namely {yfdref)o = 5.540, {y/dref) 2 b = 16.02, and {y/dref )50 = 27.115 (see figures [^. The mass ffuxes 
computed using the solution of model B are equal to rhxjdref=o = 96.28 g/s^ 'dixjdr.ef=‘ 2 b = 511.78 g/s^ and 
'dixidref=^o = 1163.77 g/s^ respectively. The different behaviour of model A reflects in the evaluation of the 
mass ffuxes which are: = 98.07 g/s, mxidr.ef= 2 b = 681.83 g/s, and rhxidref =50 = 1202.08 g/s. 

Figure shows the distribution of some relevant thermo-chemical variables at several abscissae along 
the chamber axis, comparing the results obtained by the two FPV models, an LES with reduced chemistry 
evaluation [7], a mixture-fraction based model m, and experimental results [48]. The H 2 and O 2 mass 
fraction distributions agree fairly well with LES and experimental data. The improvement with respect 
to the results of model A is more evident in the near axis region. Model A predicts a thin reacting zone 
near the burner, so that in sections at xfdref = 0.85 (see figure and xjdref = 10-8 (see figure the 
temperature distributions present a spike close to the axis. This spike characterizes also the water mass 
fraction distribution of the mixture-fraction-based model by Thibault and Boivin [49|. On the other hand, 
model B is not affected by this problem and correctly predicts the flame core. 

4 . 3 . The flame structure 

In this section the structure of the flame is studied with reference to the work of Moule et al. m and 
of Boivin et al. [7] who have analyzed the considered supersonic burner. Moule et al. m identify three 
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Figure 7: Distributions of different thermo-chemical quantities 

along the axis {y/dref = 0). Model A, red dashed line; model B, solid black line; Boivin et ah 
dashed-dotted blue line [7]; symbols, experimental data [48]. 



black line. 
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Figure 9: Chemical species distributions at xjdref = 10.8, xjdref = 21.5, and xjdref = 32.3: model A, red dashed line; 
model B, solid black line; Boivin et al. dashed-dotted blue line [?]; Thibault and Boivin dashed-dotted-dotted brown line ESI; 
symbols, experimental data m- 
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Figure 10: Model A and Model B solutions: contours of 6 jh 02 - 


characteristic regions of the flame, namely: 1) the auto-ignition zone, 10 < x/d^ef < 18, where the mixture 
prepares to ignite; 2) the stabilization region, 18 < x/dref < 26, where the flame starts; 3) the combustion 
region, 30 < xjdref < 40. 

In the first zone, the initial steps of the hydrogen oxidation take place and the reaction produces simple 
hydrogen radicals. In particular, for hydrogen combustion, the hydroperoxyl radical HO 2 can be considered 
as a good marker of autoignition. However, since the HO 2 concentration also increases in fuel-rich reaction 
zones of ignited mixtures, Boivin et al. [7] proposed a more accurate criterion for detecting autoignition. 
Such a criterion requires the simultaneous presence of high values of the normalized production rate of HO 2 
and of the reactivity of the mixture. A, defined as the positive eigenvalue of the Jacobian of the chemical 
source term associated to the autoignition chain-branching reaction [7]. For the present supersonic lift¬ 
off flame, the results of the LESs discussed in references [50] and [7] confirm the effectiveness of such a 
criterion in identifying the autoignition region as well as the following transition from the autoignition to 
the stabilization characterized by the HO 2 depletion. 

Analysing the results of the present RANS simulations, we can show that model B can provide a flame 
structure corresponding to the above scenario. Figure [Tq| shows the contours of the HO 2 production rate, 
^H 02 ^ computed using model A (left) and model B (right). This figure indicates that both models predict a 
substantial HO 2 production downstream of the jet exit, in agreement with the results of the LES by Moule et 
al. [50] (see figure 17 of their paper). Moreover, it appears that the production region is followed by depletion 
taking place in the stabilization region. Figure 11 shows the contour plot of A together with the Yho 2 ~ 
3 X 10“^ isoline, considered as the reference value for the high temperature hydrogen autoignition ISO]. It 
appears that the most reactive part of mixture is in the near-burner zone, where the mixing among the 
reactants is stronger. Both models properly evaluate the reactivity, in good agreement with the results of 
the LES of [50] (see figure 19 therein). Looking at the contours of 6 jho 2 in figure 10 and at the contours 
of A in figure El we can verify the conjecture of Boivin et al. concerning the autoignition region. In fact, 
we can see that model B predicts the autoignition at 7 < xjdref < 23 where sufficiently high values of the 
reactivity and of the HO 2 production are achieved. On the other hand, model A predicts a too fast reaction 
rate, providing a shorter autoignition region for 3 < x/d^ef <11. 

In the second zone, the flame finds its stable position and, as shown in figure it is characterized by 
high values of the heat release rate. Figure [T^ also shows the sonic isoline and the axial temperature profile 
superposed to the contours of the heat release rate. According to the LES results of m, although the 
highest temperature occurs in the subsonic region, 32 < xjdref < 38 , high values of qreact are concentrated 
in the neighbourhood of the stabilization region of the flame, 18 < xjdref < 26. For the present RANS 
computations, model A predicts high values of the temperature and of the heat release too close to the 
burner, whereas, model B provides results in very good agreement with the LES concerning both the heat 
release region and the temperature profile along the axis of the burner. 

As proposed by Boivin et al. [7], it is also possible to employ a unique parameter that identifies the 
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Figure 11: Model A and Model B solutions: reactivity contour plot with the Yho 2 isoline [3 x 10 superimposed. 




(a) Model A (b) Model B 


Figure 12: Model A and Model B solutions: heat reaction rate contour plot with the Ma = 1 isoline (black solid line) and the 
axial temperature profile (red solid line) superimposed. 
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(a) Model A 


(b) Model B 


Figure 13: Model A and Model B solutions: contour plot of a. 
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(a) Model A 


(b) Model B 


Figure 14: Model A and Model B solutions: contour plot of looh- 


occurrence of autoignition, 


a = 


(29) 


'^H02 


with positive/negative part of ujho 2 representing the production and the destruction rates oi H 02 ‘ 

According to Boivin et ah [7], the autoignition occurs when a decreases from its maximum reference value 
{oLmax ~ 0.95) to its minimum reference value {amin ~ 0.05). Figure 13 shows the two a contour plots 
obtained by model A and B, respectively, confirming that model B is capable of predicting the transition 
from the autoignition to the stabilization region at about 27 dref^ providing a good estimate of the lift-off 
height. 

shows that in the third region (30 < xjdref < 40) the combustion develops and is 


Finally, figure 
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characterized by high values of OH production. 


4.4' ^ (^'^d A distributions 

In this section, we focus on the analysis of those flow regions where la£ge differences between the pre¬ 
dictions of model A and B are observed. In particular, we compare the Psml,2 with the /d—distribution 
for Z and with the (5—function for A. All the distributions are evaluated using the values of the mean and 
variance computed by model B at three points: two of them are on the axis {y/dref = 0) and correspond 
to the lift-off heights evaluated with the two models, xjdref = 8 and xjdref = 27, respectively; the third 
point is close to the jet exit {xjdref = 0.85, yjdref = 0.8), where a spurious temperature peak is predicted 
by model A (see figures |^. 

The first point {x/dref — 8) is in the stabilization region for the solution obtained by model A and falls 
in the mixing non-burning region for the solution obtained by model B. This is clearly indicated by the 


17 





















Z = 0.34, = 0.024, A = 0.008, A"2 = 0.0 



Figure 15: PDFs of Z at the point {xjdref -,11 Idref) = (85 0) using the values obtained by model B: PsML, 2 {Z) (black solid 
line) and /3(Z) (red dashed line). 



Figure 16: PDFs of Z and A at the point {xjdrefjdref) = (27,0) using the values obtained by model B. Left panel: 
PsML, 2 {Z) (black solid line) and /5(Z) (red dashed line). Right panel: PsML, 2 {^) (black solid line) and (5(A). 
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Z = 0.21, Z"2 = 0.005, A = 0.0215, A"2 = 0.021 




Figure 17: PDFs of Z and A at the point {xjdrefIdref) = (0.85,0.8) using the values obtained by model B. Left panel: 
PsML, 2 {Z) (black solid line) and /3(Z) (red dashed line). Right panel: PsML, 2 {^) (black solid line) and (5(A). 


Table 2: Mean and variance of Z and A at three selected points. 





Z 

Z"2 

A 

7^2 

{x/dref y/dref) — 

(8,0) 

Model A 
Model B 

0.37 

0.34 

0.030 

0.024 

0.22 

0.008 

0 

{x / dref •) y / dref^ 

(27,0) 

Model A 
Model B 

0.085 

0.065 

0.0022 

0.0012 

0.56 

0.265 

0.058 

{x/dref y/dref) — 

(0.85,0.8) 

Model A 
Model B 

0.22 

0.21 

0.007 

0.005 

0.49 

0.0215 

0.021 


mean and variance values of Z and A provided in table the mean value of the progress parameter is 0.22 
for model A, whereas it is very close to zero for model B. The values of the mixture fraction are very close 
to each other. Figureshows the corresponding mixture-fraction distributions evaluated by the /d-function 
and by the Psml ,2 PDF. The two distribution of A (not shown) are both Dirac distributions since, due to 


the zero variance, the P function of model B collapses on a ^-function centred at A = 0.008. 

At the second point on the axis of the burner {x/dref = 27), the reaction is active for both models, 
as shown by the progress parameter values in table which indicates that the combustion rate is higher 
for model A. Model B also provides lower values for the mean and variance of the mixture fraction. The 
PDF of Z and A are shown in figure It appears that the Psml ,2 function is close to the /^-distribution 
due to the low variance. The main differences between the two models are clearly observed in the right 
panel of figure 16 where the Psml, 2 {^) and the S{A) are shown in order to put in evidence the dramatic 


simplification made when assuming the Dirac distribution for the progress parameter. 

Finally, considering the third point close to the burner {x/dref = 0.85, y/dref = 0.8), table shows 
that a large difference in the progress-variable mean value exists between the two models. In fact, model A 
provides A = 0.49, indicating that combustion is active, whereas the mean value of A for model B is close 
to zero. Figure 17 shows the corresponding PDF of Z and A. Again, the Psml, 2 {Z) is very close to the 
/^-distribution (left panel), whereas the Psml, 2 {A) maintains a smoother behaviour with respect to the 6- 
function adopted for model A (right panel). The latter difference can be considered a reasonable motivation 
for the smoother behaviour of the temperature predicted close to the burner by model B (see figures and 

E- 
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5. Conclusions 


This paper presents a statistical more likely distribution (SMLD) approach for the evaluation of the 
presumed probability density function (PDF) in flamelet progress variable (FPV) models for non-premixed 
combustion. The FPV model is employed in conjunction with the Reynolds averaged Navier-Stokes (RANS) 
equations. The proposed SMLD model is built evaluating the most probable joint distribution of the mixture 
fraction and of the progress variable and its adequateness and feasibility are discussed in comparison with 
the standard FPV model. The SMLD model relies on a more robust theoretical basis with a substantially 
unchanged computational cost. Although the classical formulation of the FPV approach is based on the low- 
Mach-number assumption, we solve the total energy conservation equation to improve the computation of 
compressible reacting flows. The performance of the two FPV models are discussed by analysing the results 
of the simulation of a supersonic iL 2 -Air combustion studied at the NASA Langley Research Center. This 
analysis shows that the FPV-SMLD model provides an effective improvement over the standard approach 
and is able to properly describe the flame structure in good agreement with the results obtained by highly 
resolved LES with detailed chemistry. In fact, the numerical results correctly predict the presence of the three 
characteristic regions of supersonic flames: the autoignition, the stabilization and the combustion regions. 
Moreover, the FPV-SMLD model can correctly evaluate the lift-off of the flame whereas the standard model 
is not able to predict the combustion kinetic with sufficient accuracy, providing a sudden ignition of the 
mixture close to the burner inlet. A detailed analysis of the PDF distributions at several points of the 
computational domain is provided in order to quantify and explain the differences between the two models. 
This work has shown that indeed evaluating the PDF using the SMLD approach allows one to improve the 
accuracy of the simulation of a reacting supersonic flow in the framework of RANS equations. However, 
several aspects should be considered to enhance this analysis, such as: the effects of compressibility deserve to 
be analyzed in more details employing, for instance, a compressible FPV model; the effects of the hypothesis 
that the mixture fraction and the progress variable are not statistically independent have to be studied; the 
influence of the proposed FPV-SMLD model in the LES framework should be investigated. However, all 
this topics go beyond the aim of the present paper and will be the subject of our future work. 
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